Characteristics of the wavelength of ripples on icicles 
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It is known that the wavelength of the ripples on icicles in nature is of centimeter-scale. Such 
study on morphological instability of ice- water interface during ice growth from flowing supercooled 
water film with one side being a free surface has recently been made [K. Ueno, Phys. Rev. E 68, 
021603 (2003)]. This is a first theoretical study taking into account the influence of the shape of 
the water-air surface on the growth condition of infinitesimal disturbances of the ice- water interface. 
A simpler formula to determine the wavelength of the ripples than that in the previous paper is 
derived. It seems that the wavelength of ripples is insensitive to the water supply rates, diameters of 
' the icicles and surrounding air temperatures. The details of dependence of the wavelengh of ripples 

I on these parameters are investigated. 
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I. INTRODUCTION 



Thin liquid films are ubiquitous entities in a variety of settings and display interesting dynamics. [l[ The flow of 
liquid in a thin film down a solid wall is often observed in everyday life, for example, when rain water runs in a sheet 
^Lji down a window_pane. The stability of a laminar flow of a viscous liquid running down an inclined plane was examined 
'"^ ■ by Benjamin. 2] In this case, one side of the liquid film is a free surface and the other side is a rigid plane. However, 
when the liquid flows down on a solid-liquid interface accompanying a phase transition, the solid-liquid interface may 
not remain flat if a morphological instability occurs. Such an example is ring-like ripples on icicles as shown in Fi_g^. 
[TJa), which appears when icicles are covered with a thin film of Sowing water surrounding by a cold air below 0°C. 
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^ Figure[ljb) shows the frequency of the wavelengths of ripples on natural icicles with different diameters. This is based 
. ^ on the 41 data collected from photographs which were taken in many winters at various sites in Hokkaido in Japan. 
^ ■ The growth conditions such as air temperature, water supply rate and wind speed, etc. are all different. Majorities 
^H ^l of wavelengths are within 7 ~ 10 mm and the mean wavelength is 8.5 mm. In addition to the ripple formation, as 
^ , shown by bold dashed lines in Fig. [Tfc), it is observed that many tiny air bubbles trapped in just upstream region of 
I— I ' any protruded part of an icicle move in the upward direction during ice growth. [3] This indicates that an initial fiat 
solid-liquid (ice- water) interface not only becomes unstable but also moves upwards. The similar pattern to the ripples 
on icicles can be experimentally produced on an inclined plane set in a cold room below 0°C by continuously supplying 
water from the top of the plane as schematically shown in Fig. [DJa). This experiment showed that the wavelength 
' of ripple increases with a decrease in the angle d. At 6 = tt/2, which nearly corresponds to the configuration of icicles, 
p\| ] the average wavelength of the ripples was about 7.8 mm at the water supply rate (5=160 ml/hr. 

- The fundamental building block of the morphological instability is the Mullins-Sekerka (MS) instability of a so- 
lidification front, which gives conditions for the growth of infinitesimal disturbances of a solid-liquid interface. Q 
1 According to the MS theory for a pure material, as a result of competition between destabilization due to the thermal 
T— I ' diffusion and stabilization due to the Gibbs-Thomson (GT) effect, which is the melting temperature depression due 
^-H ■ to the curvature of the interface, a pattern with a specific wavelength of Aj^f^ = 2Try/3lddo on the order of microns is 
^ ' developed. [1,0] The Aj^f^ contains two characteristic lengths. One is the thermal diffusion length = ki/V which 
is the thickness of an accumulation layer of the latent heat released by solidification. The layer is formed ahead of 
the solidification front in the semi-infinite liquid. This is usually macroscopic length determined by the values of 
the mean crystal growth velocity V and the thermal diffusivity of the liquid. The other is the capillary length 
da — TmTCpij associated with the solid-liquid interface tension F, where T,„ is the melting temperature, Cpi is 
the specific heat at constant pressure of the liquid and L is the latent heat of solidification per unit volume. This 
is a microscopic length of the order of angstroms. Later, Cristini and Lowengrub (CL) revisited the linear thoery of 
the quasi-steady diffusional evolution of growing crystals in 3-dimension, developed by MuUins and Sekerka. [1, 0] 
They reposed the problem in terms of a critical fiux rather than the critical radius in the MS theory. As a result, 
they revealed that the MS instability, that arises in the presence of constant undercooling, can be suppressed by 
maintaining fiux conditions close to critical. This near-crtical conditions can be achieved by appropriately varying 
the undercooling in time. In the MS and CL theories, the region of liquid is assumed to be semi-infinite. However, 
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since the typical thickness of flowing water film on the surface of icicles is about 100 fj,m, the thermal diffusion layer 
cannot be formed in the liquid. Furthermore, the GT effect can be neglected because the wavelength of ripples on 
icicles is about 1 cm. Therefore, we have to consider different characteristic lengths from those in the MS theory and 
to develop quite a new mechanism for the ripple formation on icicles. 
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FIG. 1: (a) Icicles with ripples, (b) Frequency of the wavelengths of ripples on natural icicles, (c) Vertical cross section of an 
iciccle with ripples. 

It is well known that the coupling of fluid flow in the liquid with phase transformation at the solid-liquid interface 
can result in the creation of new frontal instabilities, [lol ] However, the morphological instability of solidiflcation front 
in crystal growth from thin flowing liquid with one side being a free surface have never been examined. The problem 
involves both a solid-liquid interface and a liquid-air surface, in particular, both ice- water and water-air interfaces, as 
shown in Fig. EJb). No quantitative theory for predicting the wavelength of ripples has been proposed until a recent 
theoretical attempt to explain ripple formation on the surface of icicles is made. [ll|, [l3 According to their results 
of stability analysis of the solid-liquid interface, the instability of the solid-liquid interface occurs by the Laplace 
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instability due to the thermal diffusion into the air, and the instability is suppressed by the effect of liquid flow in the 
film layer. They interpreted the suppression effect as the GT like effect, which makes the temperature distribution 
in the thin water layer uniform. They concluded that ripples of centimeter-scale wavelengths appears as a result of 
competition between these two effects. However, their dependence of the wavelength of the ripples on the angle of the 
inclined plane was deviated from the experimental result by Matsuda. Furthermore, their prediction that the ripples 
move in the downward direction was opposite to the observation mentioned above. Recently, a completely different 
mechanism for the ripple formation has been proposed. [Tsl . [l^ Here we briefly summarize some of results obtained 
from our theory. 

In the first paper, [l3| the morphological instability of the solid-liquid interface of icicles was investigated by a linear 
stability analysis. We clarified that when the solid-liquid interface and the liquid-air surface are coupled even if the 
modes are sinusoidal in the linear stability analysis, it needs to take into account the influence of the shape of the 
liquid-air surface of the liquid film on the growth condition of infinitesimal disturbances of the solid-liquid interface. 
It was found that restoring forces due to gravity and surface tension, although which act on the liquid-air surface, 
is an important factor for stabilization of morphological instability of the solid-liquid interface under appropriate 
thermodyanamic boundary conditions. We determined the wavelength of ripples on the incliend plane by varying 
the angle 9 at the same water supply rate per width as the experiment by Matsuda. Our theoretical result showed 
that the wavelength increases with a decrease in the angle, which was in good agreement with experimental result 
(see FIG.4 in Ref. Tsl). Moreover, our theory predicted that the solid-liquid interface moves in the upward direction, 
which is consistent with the observation in Fig. [TJc). 





FIG. 2: (a) Schematic view of ice growth from a thin layer of supercooled water flowing down an inclined plane to the horizontal 
by angle 9. Q ml/hr is the water supply rate from the top and / cm is the width of the gutter. Shaded portions within the 
gutter show nearly periodic ripples of ice. (b) Schematic view of vertical cross section {x, y) of the inclined plane. The thick 
solid line is the ice-water interface and the soUd line is the water-air surface, ho and uq are the mean thickness and the surface 
velocity of flowing liquid film in unperturbed state. 

In the second paper, [iJl the magnitude of restoring forces due to gravity and surface tension changes depending 
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on the wavenumber of the disturbed soHd-hquid interface. As a result, the shape of the Uquid-air surface is changed. 
There is a phase difference (j) between a disturbed sohd-hquid interface and a distribution of heat flux at the solid- 
liquid interface, whose magnitude is the cause for destabilization and/or stabilization of the solid-liquid interface. We 
showed that if — 7r/2 < </> < 0, a disturbed solid-liquid interface is unstable to form ripples, while if — tt < cj) < —tt/2, 
a disturbed interface is stable and becomes flat. We also showed that the sign of (j) is related to the direction of 
movement of the solid-liquid interface. Since < 0, the solid-liquid interface moves in the upward direction. 

Our dispersion relation contains the mean thickness of the liquid film, ho = [2i^/{gsme)]^/^[3Q/{2l)]'^/^, where Q 
is the water supply rate, I is the span width of the gutter (see Fig. [H^a)), v is the kinematic viscosity of the liquid and 
g is the gravitational acceleration. The experiment by Matsuda was performed for various temperatures of the cold 
room and the water supply rates under fixed span width / = 3 cm. After trial and error, he found that ripples were 
clearly formed at the room temperature of —8 ± 1°C and the supply rate Q — 160 ml/hr. [3] Hence, in the previous 
papers, IT^ we determined the wavelength of the ripples by varying the angle 9 at the same water supply rate 
per width, Q/l = 160/3 [(ml/hr)/cm], as the experiment by Matsuda. However, the equation of ho above indicates 
that the thickness also changes by varying the value of Q/l- Therefore we can expect the wavelength of the ripples to 
depend on Q/l. Nevertheless, it was claimed that the wavelength of the ripples on natural icicles is independent of 
water supply rate and diameter of the icicles because the amount of Q is expected to vary with the radius R of icicles 
to keep the ratio Q/l = Q/{2TrR) during icicle growth, [llj In order to clarify it, we investigate the dependence of 
the wavelength on Q/l and ambient temperature. We will present a simple formula to determine the wavelength of 
ripples. We also predict the dependence of the wavelength of ripples on the surface tension of the liquid-air and the 
dependence of the movement velocity of the solid-liquid interface on Q/l. 



II. THEORETICAL FRAMEWORK 

We give the theoretical framework developed in the previous papers p^fl^ with some modification and supplement. 
As schematically shown in Fig. [2{a), the water at above nearly 0°C is supplied from the top nozzle at the rate Q 
ml/hr and flowing down on an inclined plane with span width I cm. Q Since the apparatus is set in a cold room 
below 0° C, the water becomes in the supercooled state as it flows down the plane. Ice grows from a portion of the 
supercooled water layer through which the latent heat of solidification is released into the ambient air below 0°C. 
The rest of the water drips from the lower edge of the plane. Our analysis is restricted to two dimensions (x, y) of 
the vertical plane, and we assume that the region of the crystal is semi-infinite as shown in Fig. [IJb) because the 
thickness of the liquid film is very thin compared to the typical redii of icicles. The x axis is parallel to the inclined 
plane and the y axis is normal to it. The position of y = is the solid-liquid interface in the frame of reference moving 
at the solid-liquid interface with an unperturbed solidification velocity V . 

The velocity profile in the liquid film in the unperturbed state is given by the parabolic shear fiow: [2, [isj 



C7(,) = uo 2f- f , (1) 




which is parallel to the x axis as shown in Fig. EJb), where ho is the mean thickness of liquid film and uo is the surface 
velocity of liquid film in the unperturbed state. Then the water supply rate per width, Q/l, is given by 

Q/l = Uo (^2^ - fj) = 

Since uo = gh"^ sin 9 / {2v), Ts'] ho and uq are not independent quantities. Therefore, it is convenient to express them 
by experimentally controllable parameters Q/l and 9: 

The following material properties of water are used: the kinematic viscosity, ly = 1.8 x 10^^ in^/s; the density, 
pi = 1.0 X 10'^ kg/m'^; the thermal diffusivity, Ki = 1.3 x 10~^ m^/s and the specific heat at constant pressure, 
Cpi ~ 4.2 X 10^ J/(m^°C). The values of the surface tension of the water-air surface, the latent heat of solidification 
per unit volume and the ratio of the thermal conducutivity of ice to the thermal conducutivity Ki of water are 
7 = 7.6 X 10~^ N/m, L = 3.3 x 10® J/m'^ and n = Ks/Ki=3.92, respectively. Artificial icicles were made in a 
cold laboratory room by controlling temperatures, water supply rates and wind speeds. [3] The water supply rates 
were varied in the range of Q = 14.4 ~ 108 ml/hr. Then the corresponding dripping rates from the tip of icicles 
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were 15 ~ 2 seconds. On the other hand, the water supply rate in the experiment by Matsuda [j] was kept at only 
Q = 160 ml/hr. In the latter case, the values of /iq and Uq si 9 — 7r/2 and / = 3 cm are 93 /im and 2.4 cm/s, 
respectively. Figure [3] shows the dependence oi on Q/l ior = tt/2 and n/18. ho increases only gradually with 
Q/L As dimensionless numbers, we define the Reynolds number. Re = uoh^/v = / {2v)]{Q / 1); the Peclet number, 
Pe = uoho/Ki = [3/{2k,i)]{Q/1); the Froude number, F = uo/ighoY^^ = {sin6l/(2i^)}i/2{3g/(2Z)}i/2 and the Weber 
number, W = j/{pihoul) = j /[pi{g sinO / {2v)y^^{3Q / {2l)Y^^]. We note that Re and Pe depend on only Q/l, while 
F and W depend on both Q/l and 6. For the same value Q/l =160/3 [(ml/hr)/cm] as used in the experiment by 
Matsuda and 9 = 7r/2, these dimensionless numbers take the following values. Re ~ 1.23, Pe ~ 17.1, F ^ 0.79 and 
W ~ 1.42 X 10^. 
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FIG. 3: Dependence of the thickness ho of the liquid film on the rate of volumetic flow per width, Q/l [(ml/hr)/cm] for 8 = 7r/2 
and 7r/18. 



A. Governing equations 



The equations for the temperature in the liquid T;, solid and air Ta are given by 

dTi dTi dTi / d^Ti &^ 
dt ^ dx ^ dy '^^ \ dx'^ dy-^- ' 



dTs - dTs _ /a^T, d^Ts 



= 1^ + -TOT ' (6) 



dt dy \ dx^ dy 

where t is time, k;, and are the thermal diffusivities of the liquid, solid and air, respectively. In Ref. [Tsl , u 
and V were the velocity components in the x and y direction measured in the laboratory frame, here we replace them 
with the velocity components measured in the frame of reference moving at the solid-liquid interface with V. An 
experimental observation shows that the increase in the wind speed in the air leads to the increase in V. 3 We will 
show that the change in V affects the magnitude of amplification rate for disturbances of the solid-liquid interface, 
but the characteristic wavelength of the ripples does not change. As far as we are concerned with determining the 
wavelength, it is sufficient to assume no wind in the air. 
The Navier-Stokes equations for u and v are 

du du du 1 dp / d^u d^u\ . ^ 
dt dx dy pi dx \ dx"^ dy'^ ) ^ 

dv ^ dv ^ dv I dp ^ / d'^v ^ d'^v\ ^ 
dt dx dy pi dy \dx^ dy"^ J 
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where p is the pressure, pi is the Uquid density. The equation of continuity is 

du dv ^ 
dx dy 

From (j9]), u and v can be expressed as u = uodip/dy and v = —uodip/dx, where tp is the stream function. 



(9) 



B. Boundary conditions 

1. Thermodynamic boundary conditions 
The continuity of the temperature at a disturbed solid-liquid interface y = (^{t, x) is 

The heat conservation at the solid-liquid interface is 

dTi 



y=C oy 



(10) 



(11) 



where L is the latent heat per unit volume, and and Ki are the thermal conductivities of the solid and liquid, 
respectively. On the other hand, the continuity of the temperature at a disturbed liquid-air surface y = ^(t,x) is 



(12) 



where Tia is a temperature at the liquid-air surface, which will be determined below. The heat conservation at the 
liquid-air surface is 



dy 



-K, 



y=i 



dy 



y=i 



(13) 



where Ka is the thermal conductivity of the air. We refer to these as the thermodynamic boundary conditions. 



2. Hydrodynamic boundary conditions 

We neglect the density difference between ice and water from the first. Then both velocity components at the 
solid-liquid interface must satisfy 



The kinematic condition at the liquid-air surface is 
At the liquid-air surface the shear stress must vanish: 



du 
dy 



dv 

y=i dx 



y=i 



= 0, 



(14) 



(15) 



(16) 



and the normal stress including the stress induced by surface tension 7 of the liquid- air surface must just balance the 
pressure Pq of the atmosphere: 



y=C ' dx 



(17) 



We refer to these as the hydrodynamic boundary conditions. 
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C. Perturbations 



As seen in Fig. [Ija), ring-like structure encircles the icicles and there is no noticeable azimuthal variation on the 
surface of the icicles. Therefore, it is sufficient to consider only a one dimensional perturbation in the x direction of 
the solid-liquid interface: 



C,{t, x) — Cfe exp[o'< + ikx], 



(18) 



where k is the wavenumber and cr = (Xr + ^f^, with Ur being the amplification rate and Vp = —Ui/k being the phase 
velocity of the perturbation, and Cfe is a small amplitude of the solid-liquid interface with the wavenumber k. We 
separate ^, T/, T,, T^,, tp and p into unperturbed steady fields and perturbed fields with prime as follows: C = + d 



Ti=Ti+ T' T, 



T, 



T' 



T 



ijj = ip + Tp' and p = P + p'. We suppose that the perturbations of the 



liquid-air surface, the temperature in the liquid, solid and air, the stream function, and pressure are expressed in the 
following forms: 



Tlit, X, y) 

Ta{t,x,y) 
■ip'{t,x,y) 
\ P'{t,x,y) ) 



( \ 

giiv) 
9s{y) 
9a{y) 

F{y) 
V n(y) J 



exp[cri -|- ikx], 



(19) 



where ^fe, gi, g^, ga, F and 11 are the amplitudes of respective perturbations and they are assumed to be of the order 
of Cfc- 

The following calculations are based on a linear stability analysis taking into account only the first order of C,k- 
Furthermore, we make two approximations. First, we use the long wavelength approximation, (T], [^j which is valid 
when the fluid film thickness is much less than the characteristic scale of ripples. Then we define a dimensionless 
wavenumber /i = kh^ and neglect higher orders in ^. Second, we use the quasi-stationary approximation. 0] We 
can neglect the time dependence of perturbed part of temperature in each phase and that of liquid flow because these 
fields respond relatively rapidly to the slow development of the perturbation of solid-liquid interface. 



D. Equation for the amplitude of perturbed part of the stream function and its solution 



When we substitute ip' and p' into the equation for the perturbed parts of ([7]) and ^ and 11 is eliminated from them 
by cross differentiation, we obtain the Orr-Sommerfeld equation for the amplitude F{y^) = f(jJ*)Ck of the perturbed 
part of the stream function: 



dyt 



ifiRe U2y, - yt)-^ + 2f 



dyl 



(20) 



where y^ — y/fiQ. Since the values of /ii?e is very small, we neglect the term on the right hand side of ([20)1 . Using the 
following linearized forms of hydrodynamic boundary conditions (|14p - (jl7p . 



dy* 



f\v, = lCk 



dyl 



/l..=o=0, 
Cfc = 2^fc, 
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dyl 



Cfc - iOL^k = 0, 



the solution of d^f/dyf = is given by 



f{y*) = -'^y* + 



3(2 — ia) 2 



6 — ia 



6 — ia 



(21) 
(22) 
(23) 

(24) 
(25) 
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FIG. 4: (a) Dependence of the amplitude ratio |^fc]/Cfc on the wavenumber k for 9 = 7r/2 and n/lS at Q// =160/3 [(ml/hr)/cm]. 
A typical value of the capillary length a of water is about 2.8 mm, which corresponds to the value k = 2244 indicated by 
the arrow. The shape of the water-air surface (solid lines) for a disturbed ice- water interface (thick solid lines) with a amplitude 
Ck/ho = 0.1 when the wavelength (b) 1 cm (k ~ 628 m~^) and (c) 2 mm (k ~ 3142 m~^) for 9 = it/2, (d) 1 cm and (e) 2 mm 
for 6 — 7r/18. In the case of (b) and (c), ho ~ 93 ^m, while for (d) and (e), ho « 168 pm. 



where 



, ^i^ReW ^2 cot dhok-\ a'^h^k^ (26) 

sin 



represents the restoring forces due to gravity and surface tension [2| and a = [y /(pig)]^^^ is the capihary length 
associated with the surface tension 7 of the liquid-air surface. [15| Since the value of W is very large as estimated 
above, the /jP term in a cannot be neglected. 

It should be noted that the difference between the stability analysis of a laminar flow and that of a solid-liquid 
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interface. In the stability analysis of the laminar flow of a viscous liquid flowing down an inclined rigid plane, |2'| the 
boundary conditions for the flat plane, u|j,=o = and w|y=o — 0, are used instead of p^ . These and ([T5]) and 
are used to determine four constants arising from the solution of the fourth order differential equation ()20|) without 
neglecting the term of time derivative of /. The last boundary condition ([T7| then be used to find the dispersion 
relation for the disturbance of the free surface. In the present case, however, there are five boundary conditions in 
([T4| - (|17p . From (fTS]). the relation between the amplitude of the liquid-air surface and the amplitude Cfc of the 
solid-liquid interface is obtained as follows 

= -f\y,=iCk = TT^Cfc. (27) 

Figure m^a) shows the dependence of the amplitude ratio |$fc|/Cfe on the wavenumber fc at Qjl ~ 160/3 [(ml/hr)/cm] 
for Q = Tr/2 and ir/lS. |^fc|/Cfe decreases monotonically with k as shown in Fig. IH^a). 

In the low wavenumber region in Fig. Hl^a) , the restoring forces due to gravity and surface tension on the liquid-air 
surface can be neglected because the value of a is very small, therefore, the amplitude of the liquid-air surface is 
almost the same as that of the solid-liquid interface as shown in Figs. lUJb) and (d). In this case, fluid flows down 
along the disturbed solid-liquid interface. On the other hand, in the high wavenumber region in Fig. SJa), since the 
value of a is large, the shape of the liquid-air surface tends to be flat against the solid-liquid interface as shown in 
Figs. Sfc) and (e). We can see that fluid flows down without feeling the shape of the solid-liquid interface. It should 



be noted that in the Ogawa and Furukawa's theory llj the amplitude of the liquid-air surface is always the same as 
that of the solid-liquid interface irrespective of the wavenumber of the solid-liquid interface because the effect of the 
restoring forces on the liquid-air surface is not taken into account. 



E. Equation for the amplitude of perturbed part of the temperature in the liquid and its solution 

Solving the equation for the unperturbed part of @ with the boundary conditions T/ = T„i at y = and T/ ~ Tia 
at y = /iQ, the solution is Ti{y) — T„i — Giy, where Gi = —dTi/dy\y^o = {Tm — Tia)/hQ is the unperturbed part of 
temperature gradient in the liquid fllm and G; > means that water is in supercooled state. When we substitute the 
assumed forms of T/ and ip' into the equation for the perturbed part of (|3|) , we obtain the equation for gi : 



-j-^ - i^J.Pe{l - z^)gi = inPef{z)GiCk, 

where we have introduced a new varaible, z = \ — y/h^. 

When putting the right hand side of (j28p equal to zero, the independent solutions are given by 



(28) 



bi{z) = exp ( -i(-i/iFe 



,1/2^2 



1^1 (J 



1 





1 








^4 



(-^/iFe)l/2 
ifiPe 



(-^^Fe)l/2 



(29) 



(30) 



where iFi is the confluent hypergeometric function. 16] Noting that the Wronskian of the two solutions 4>i{z) and 
<j>2{z) equals to unity, the general solution of (|28l) can be written as 



gi{z) = BiM^) + B^M^) + if^PeGiCk / {Mz)M^') " M^)M^')} f{z')dz\ 

Jo 



(31) 



where Bi and B2 are unknown constants to be determined from the thermodynamic boundary conditions at the 
liquid-air surface. 

With the boundary conditions Ta — Tia at y — ho and Ta — Too s-t y — ho + S, the solution of the equation for the 
unperturbed part of ([6]) is 



Taiy) = Tia - Ga{y - kg) , 



(32) 



where Ga = —dTa/ dy\y=ho = iTia — T^)/6 is the unperturbed part of the temperature gradient in the air at the liquid- 
air surface, and we regard 6 as the thickness of a thermal boundary layer. Strictly speaking, the thermal boundary 
layer, which is created around icicles warmer than surrounding air, can be provided from similarity solutions for the 
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coupled Navier-Stokes and heat transport equations in the Boussinesq approximation. However, we assume that the 
solution (1321) is different from the exact form of the similarity solution only by the multiplication of an order one 
constant. |l7| Therefore, we estimate the value of 5 from Eq. (4) in Ref. 

Substituting the assumed form of into the equation for the perturbed part of ([5]), we obtain 



9a{y) = Tkaexp[-k{y - ha)] 



(33) 



where Tka is the amplitude of perturbed part of the temperature in the air. Linearizing (|12p at y = ho gives to the 
zeroth order in 



Tl\y—ho — Ta\y=ho — Tia, 

and to the first order in ^fc, noting that gi\z=o — Bi since 4>i\z=o — 1 and (l)2\z=o = 0, 

Linearizing (IT^ aX y — ho gives to the zeroth order in 

Kid = KaGa, 

and to the first order in ^k, noting that dgi/dz\z=o = B2 since d(f>i/dz\z=o — and d(f>2/dz\z=o = 1, 

K1B2 = flKaTka- 

From ([55)) -(|37 p . Bi and B2 are determined and ([5T|) becomes 



-f\z=o{Mz) + ^^Mz)} + ^^^Pe / {Mz)Mz') ~ Mz)Mz')} .f{z')dz' 



where we have used the relation — — /|z=oC/c- Hence, (|28p can be written as 

d^Hiiz) 



GiCk = Hi{z)GiCk, 



dz^ 



i^iPe{l - z^)Hi{z) = ifiPef{z). 



(34) 
(35) 
(36) 
(37) 

(38) 
(39) 



It is found that the perturbed part of heat transport in the liquid film in the long wavelength region is dominated by 
mean shear flow U of the second term on the left hand side and normal flow v' generated from a disturbance of the 
solid-liquid interface on the right hand side. The heat transport by thermal diffusion can be neglected because this 
term appears from the second order in /i. 



F. Dispersion relation 



Linearizing (jlip at y — gives to the first order in Cfe, we obtain 



LaQk = Kg—— 
dy 



v=o dy 



(40) 



With the condition that the disturbance in the solid must vanish far from the solid-liquid interface, the solution of 
the equation for the perturbed part of ^ is given by 



9s{y) = Tksexp{ky), 



(41) 



where Tks is the amplitude of perturbed part of the temperature of the solid, which is determined from the perturbed 
part of pU]) as follows 



GlCk + gi\y=0 — Tks 



(42) 



Substituting ([551) and ([1T|) into (HO]), the dispersion relation for the disturbance of the solid-liquid interface is 
obtained: 



V {d_m 

dz 



(43) 
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where n = Ks/Ki, and 



Jo 



and 



dHi 
dz 



=1 dz 



We can express and (|30p by the expansion in terms of fiPe as follows 

Mz)^l + i( -z" - —z^\ iiPe + {-—z" + —z^ - — {fxPef + 
' \2 12 V 24 360 672 / ' 



(44) 



(45) 



(46) 



(j)2(z) =z + i(-z^- —zA aPe + ( —z"^ + —z'^\ iuPef + • • • . (47) 

' V6 20 y ^ V 120 2520 1440 ) > ^ ^ ' 

Putting z = 1 in dUl) and dH]), we obtain = 1 + hi{^JLPe)/l2 - 239(Ai-Pe) ^10080 + ••• and 02^=1 = 

1 + 7i{fiPe)/60 — 13(/ii-'e)^/3360 + The ratio of the term of second order in fxPe to the first one in 

(pilz^i and 4'2\z=i are about 5.7 x lO^^^Pe and 3.3 x 10~^^Pe, respectively. When the value of /iPe = 
[34/3/(2K;)][j//(5sin6')]i/3(Q/;)4/3(27r/A) takes about 100, the first term and second one become the same order. 
This is possible when Q/l ~ 1000 [(ml/hr)/cm] for typical wavelength A of ripples on icicles. For natural icicles with 
radii R, Q/l does not necessarily take a definite value. Since the typical values of Q over icicles are on the order of 
tens of ml/hr and their radii are usually in the range of 1 ^ 10 cm, Q we vary the value of Q/l over the range of 
10 ^ 100 [(ml/hr)/cm]. As far as we are concerned in such a range of Q/l, fJ,Pe ~ 0(1) for the wavelength of ripples 
on icicles, then it is sufficient to consider up to the first term in /xFe of (|46| and (|47)) in the following calculations. 



III. RESULTS 



The amplification rate ct^ and the phase velocity Vp — —ai/k are obtained from the real and imaginary part of (j43p . 
respectively, as follows: 



V 
ho 



vi^iPe) + {36 - §a(/iPe)} ~^a{fiPe) - a'^ + fx {36 ~ ^a{fiPe)} 



36 + a2 



36 + a2 



(48) 



V 



a2(/iPe) + /i{6a + 9(AtPe)} 6a- ^a^i^iPe) + |I{6a+'^{^J^Pe)} 



With the assumption of = T„i in the solid, linearizing (jlip at y = gives to the zeroth order in C,k, 

CplKl{T„i — Tia) 



(49) 



V = 



Lha 



(50) 



where Cpi is the specific heat at constant pressure of the liquid. From ([50|. the limit F/iq/k; ^ 1 corresponds to 

Cpl {Tm — Ti 



< 1. 



(51) 



This indicates that the supercooling Tm — Tia in the liquid is much smaller than L/Cpi. In this case, thermal diffusion 
layer of the thickness ki/V ahead of the solidification front in the liquid is not formed because ^ ni/V, and the 
latent heat is released into the air through the thin liquid film. For the typical values V ~ 10^^ m/s and ho w 10^** 
m, this condition is indeed satisfied. 
From dSH), Tia in (jH]) is obtained as 
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FIG. 5: (a) The dependence of V on Q/l for various air temperatures Too. S is the thickness of the boundary layer around 
icicles, (b) Amplification rate ar vs the wavenumber k at Q/l = 50 [(ml/hr)/cm] and Q = -k jl for various Too- 



If the values of Too, and 8 are given, the value of Tia is determined. Substituting ([52|) into ([50| yields 



if, i5 • 

ifa ha 



Lho 1 

Since Ki/Ka « 22.7 » 1 and (5 » /iq, ([5^ can be approximated as 

?(' T — T 



V: 



L 



(53) 



(54) 



Indeed, we estimate the value of V from (|53l) for various Q/l and Too using Eq. (4) in Ref. [iTl when C = 1 and z = 0.1 
m. The value of V increases with a decrease in the ambient air temperature Too as shown in Fig. [DJa). However, it 
is found that V remains constant against Q/l, hence V does not depend on ho as indicated in ([54]). This behavior 
agrees with an observation that growth velocity of diameter of icicles is almost constant with an increase in water 
supply rate (see Fig. 6b in Ref. |3|). 

Figure [5jb) shows the amplification rate ar for disturbances of the wavenumber k of the solid-liquid interface. It is 
found that the change in V by ambient air temperature affects the magnitude of maximum point of ar, but does not 
make change the characteristic wavelength Amax of the ripples, which is defined at the wavenumber at which ar takes 
a maximum value. Heat transport through the air surrounding icicles can be greatly influenced by the presence of 
forced convection. Such an enhancement of heat transport around growing icicles increases the mean growth rate V 
and results in early appearence of ripples because the value of l/cr becomes shorter without changing the wavelength 
of ripples. This theoretical prediction should be checked by experiment of icicle growth under both calm environment 
and forced convection. 

Figure[6l(a) shows the dependence of Amax on Q/l in the range of 10 ^ 100 [(ml/hr)/cm] a.t 6 ^ tt/2. The thick 
solid line is the result obtained from (|^5| . Since fiPe ^ 1 and a ~ 1 for the wavelength of ripples on icicles, (US]) can 
be approximated as 



V_ 

ho 



36" 



36/i 



(55) 
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FIG. 6: (a) Dependence of Amax on Qjl [(ml/hr)/cm] at = 7r/2. Thick solid line: Amax obtained from (|48|l . thick dashed line: 
Amax ~ 0.17 X (Q/Z)'*''^ (cni) obtained from the formula Amax = 2-K(cP'h,aPel'Sf'^'^ . Two solid circles are experimental results 
by Matsuda. [J] (b) Dependence of Amax on the surface tension 7 at Qjl = 50 [(ml/hr)/cm] and Q — 7r/2. Thick solid line: 
Amax obtained from l|48p . thick dashed line: Amax ~ 2.31 x 7^^^^ (cm) obtained from the formula Amax ~ 2-K{a^hoPe/?>)^^^ . (c) 
Dependence of Vp on Q/l aX 6 = 7r/2 when V = 10~^ m/s. 



Here we have neglected the second term in (|48p because this term is smaller than the first one and the wavenumber 
at which ar takes a maximum value is almost the same without the second term. However, we note that this formula 
holds near 9 = tt/2 because the value of a increases with a decrease in 9. From dar/dk = in ((55|) . we obtain a 
simple formula to determine the wavelength of ripples: 



1/3 



0.17 



4/9 



(cm). 



(56) 



14 



This result is shown by the thick dashed Une in Fig. IUJa), which gives a good approximation to the thick sohd 
hne. Amax increases only gradually with Q/l. Indeed, Amax = 5 13 mm for Q/l = 10 ~ 100 ml/hr. These 
wavelengths are within the distribution of the wavelength of ripples on natural icicles in Fig. [Ijb). The numbers 
17 min, 28 min and 38 min attached near the thick solid line represent the characteristic time of appearence of 
ripples for Q/Z=10, 50 and 100 [(ml/hr)/cm], respectively. These values have been estimated from the value of 
1/CTr at the most unstable mode when V ~ 10^® m/s. It is found that the characteristic time increases with an 
increase in Q/l. On the other hand, in the Ogawa and Furukawa's theory the formula to get the wavelength of 
ripples is given by A^ax ~ {2t: /2.2)hoPe w 2.3 x 10^3(g/^)4/3 -y^ ^^As leads to A^ax = 0.5 - 11 mm for 

Q/l = 10 ^ 100 ml/hr. However, the two formulae are obviously different. Each Amax shows completely different 
power law behavior with respect to Q/l. A,„ax obtained from our theory contains two characteristic lengths, the 
capillary length a — [y / [pig)]^^"^ and the mean thickness = [2v/{gsm.e)Y/^[iQ/{2l)Y'^ of the liquid film, and 
one dimensionless number Pe = [3/(2k/)]((5/^). On the other hand, Amax obtained from the Ogawa and Furukawa's 
theory contains only one characteristic length /ig and one dimensionless number Pe. Two solid circles in Fig. Eta) 
are experimental result by Matsuda. [3| One is the mean wavelength of 7.4 mm for ripples formed on a wooden round 
stick with the diameter of 1.5 cm, and the other is the mean wavelength of 8.3 mm for ripples formed on a vertical 
plane within the gutter with the width 3 cm. In both experiments, water supply rate from the top was fixed at 160 
ml/hr. There is no experimental data for other Q/l. Matsuda mentions that ice itself does notgrow when the water 
supply rate is too small or too large. [4i] The similar thing is also mentioned by Makkonen. [18| These suggest that 
there may be an appropriate range of Q/l to clearly form ripples. The present theory cannot explain this fact because 
the growth rate V is constant even in the case of too small and too large water supply rate as shown in Fig. [SJa). In 
order to confirm our theoretical prediction for the dependence of Amax on Q/l, it needs to do experiments to see the 
behavior of the growth rate of ice for much water supply rate than that shown in Fig. 6b in Ref. |3| and to modify the 
present theory which takes into account the dependence of the mean growth rate of the radius of icicles on the water 
supply rate. 

In our theory, the capillary length associated with the surface tension of the liquid-air surface is one of the im- 
portant characteristic lengths in determining the wavelength of ripples. Our theory predicts that the wavelength of 
ripples becomes shorter with the reduction of the surface tension as shown in Fig. IH^b). The thick solid line is the 
result obtained from (l48l) . The thick dashed line obtained from the formula Amax — 2n{a^hoPe/3y^'^ gives a good 
approximation to the thick solid line. It needs to measure the wavelength of ripples when reducing the surface tension 
of water by adding a surfactant. 

Figureinic) shows the dependence of the phase velocity \vp\ on Q/l at 6 — tt/2 and V — 10^^ m/s. Here we define 
the value of Vp at the wavenumber at the maximum point of ar. At this wavenumver, (|49|) takes a negative value. 
It shows that the change in Vp for various Q/l is small, which takes almost the same order as the mean growth rate 
V = 10~^ m/s=0.36 cm/hr. The negative sign means that the solid-liquid interface moves against the primary shear 
flow, which is consistent with the observation that many tiny air bubbles are trapped in just upstream region of any 
protruded part of the icicle and line up in the upward direction during icicle growth as shown by bold dashed lines 
in Fig. [IJc). [3*1 It needs to measure precisely the magnitude of velocity of the motion of the ripples under carefully 
control of Q/l and air temperature. 

IV. SUMMARY AND DISCUSSION 

We have investigated the morphological instability of the ice- water interface during ice growth from flowing super- 
cooled water film with one side being a free surface. We have obtained a simple formula, Amax — 27r(a^/io-Pe/3)^/'^, 
to determine the wavelength of ripples. Amax depends on the two characteristic lengths, the thickness ho of the water 
film and the capillary length a associated with the surface tension of the water-air surface, and one dimensionless 
number Pe, which is the ratio of heat transfer of liquid fiow to that of thermal diffusion. It should be stressed that 
Amax includes the fiuid effect in the thin water film througth the parameter Q/l. This is contrast to Amfx = 27rV^dC?o 
obtained from the MS theory, where Id is the thermal diffusion length and do is the capillary length asoociated with 
the solid-liquid interface tension. It is found that each characteristic length scale is completely different. Icicles in 
nature possess a wide range of Q/l. Nevertheless, the wavelength of the ripples on icicles seems to be insensitive to 
the supply rate Q and the radius R of icicles. This is not because Q varies with R to keep the ratio Q/R as discussed 
in Ref. Ull, but because Amax does not strongly change for the range of typical values of Q/R on icicles. Indeed, we 
have found that Amax increases only gradually with Q/l. 
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